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Abstract 

The phase diagram of quantum electron bilayers in zero magnetic field is ob- 
tained using density functional theory. For large electron densities the system 
is in the liquid phase, while for smaller densities the liquid may freeze (Wigner 
crystallization) into four different crystalline phases; the lattice symmetry and 
the critical density depend on the the inter-layer distance. The phase bound- 
aries between different Wigner crystals consist of both first and second order 
transitions, depending on the phases involved, and join the freezing curve at 

three different triple points. 
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Coupled electron layers, such as those realized in semiconductor double quantum wells, 
are a system of current interest exhibiting a lot of new physics not present in single layers 
(SLs). The focus is on the peculiar properties stemming from the Coulomb interaction 
between the layers, such as the recently reported even denominator fractional quantum Hall 
states In this paper we are concerned with the influence of the inter-layer interaction on 
Wigner crystallization in electron bilayers (BLs) in zero magnetic field. 

At sufficiently low densities a degenerate liquid of charges freezes into a lattice, which 
is called a Wigner crystal (WC), as the kinetic energy cost of the localization is more than 
compensated by the gain in Coulomb energy. Wigner crystallization has been investigated 
intensively over the past few years in SLs reahzed in semiconductor heterostructures; because 
a WC forms at low density, disorder-induced localization has been a serious competing 
mechanism which has prevented an unambiguous identification of this state. Normally, high 
magnetic fields are used to quench the kinetic energy and favour the crystallization, and a 
SL WC is believed to form for filling factors < 1/5 0. 

It has been suggested that in bilayer structures localization occurs at higher densities 
than in SLs, making this an appealing system for the investigation of the WC. The reason 
is that each layer acts as a polarizable background for the other, favouring the formation of 
an inhomogeneous phase |^ . Evidence of the formation of a BL WC in high magnetic fields 
has been found recently Q at larger filling factors (z/ < 1/4 per layer) than in SLs. On the 
other hand, in contrast to the SL WC, where the triangular lattice has the lowest energy 0, 
the BL WC might have several different lattice structures, as a result of the competition 
between inter- and intra-layer interactions, leading to a complex phase diagram A 
minimization of the classical potential energy |^ shows that five different crystal lattices are 
favoured in different ranges of inter-layer distance d and charge density n. In the classical 
regime the static energy and, therefore, the stable structure depends on the dimensionless 
product d^/n] in contrast, in the quantum regime the situation is more complicated as, due 
to the kinetic energy contribution, d and n do not scale out. 

A comprehensive picture of the phase diagram of the BL WC is still lacking. The 
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quantum Hall regime was investigated recently within the Hartree |^ and the Hartree- 
Fock approximation ||^ as a function of the filling factor, the layer separation, and the 
inter-layer tunneling. In the quantum Hall regime several solid phases have been identified 
which are analoguous to the classical ones [0], as it is expected from the fact that the 
high magnetic field quenches the kinetic energy and leads the system towards the classical 
regime. However, the question whether all phases survive quantum fluctuations could not 
be addressed. Indications of the freezing transition in zero magnetic field was found within a 
linear response theory P] , but no information on the stable lattice structure in the different 
ranges of d and n was obtained. 

In this paper we calculate the phase diagram of two coupled electron layers in zero 
magnetic field by density functional theory (DFT) ||10[. The successful application of DFT 
to the freezing transition in several systems makes this method appealing, in view of its 
formal simplicity and the possibility to include accurately [within the local density approx- 
imation (LDA) |Ty]] exchange and correlation contributions to the energy, once these have 
been calculated (e.g., from simulations) for the liquid phase. DFT allows us to study both 
homogeneous and inhomogenous phases on equal footing, and, therefore, to identify both 
solid/liquid and solid/solid phase boundaries. The system is modelled by two layers of 
mobile electrons with the same average density [0, compensated by a fixed uniform back- 



ground of positive charge. Our results can be summarized as follows: we find that four 
different crystalline phases can freeze from the liquid, depending on d and n. Below the 
freezing density, structural transitions between the different lattices occur, the transition 
being of the first or second order, depending on the phases involved. Three triple points 
separate the liquid and the solid phases. 

To set up a DFT, one has first to identify a functional of the one-particle charge density 
p(r) which is a good approximation of the exact free energy of the many-body system of 
interest which, of course, is in general unknown |jl2[. Recently, a functional was proposed 
for electrons in a SL ||T3[ which, in terms of the plasma parameter = {y/wnaB)^^ {clb is 
the Bohr radius) predicts freezing at Vg ~ 32, in good agreement with the value = 37 ± 5 



obtained from simulations |T^. This approach [0 proceeds by writing the density matrix 



p(r, r') in terms of the single-particle density p(r) with the ansatz 

p(r,r') = /(r,r')[p(r)p(r')]^/'. (1) 

The function /(r,r') can be calculated exactly for the non-interacting electron liquid, 
/(r,r') = 2Ji{kF\r — r'D/kplr — r'\, where Ji(x) is the Bessel function of the first kind 
and kp = (27rpo)^^^ is the 2D Fermi wavevector, with po the electron liquid density. Then 
one writes the kinetic energy , and the exchange energy (units of energy are Rydbergs) 

Ek = -J cir VrVr'p(r,r')|r=r', = -(1/2) JdrJ rfr'|r - r'^Xr, r')|^ (2) 

in the LDA by defining a local Fermi wavevector kp^v) = (27rp(r))^/^. Inserting (|l|) in Ek 



and Ex one obtains the total energy of a SL |jT5| with a (in general inhomogeneous) charge 
density p(r) 

Eo[p] =njdr [p{r)f + (1/4) J rfrp(r)-i| Vrp(r)p + J dr' J dr\r - [p(r) - Po] [p(r') - Po] 
-m)^ I dr [p{r)f' + / cirp(r)e,(p(r)). (3) 

The first two terms stem from E^, the third term is the direct intra-layer Coulomb inter- 
action, with Po the average charge density, and the last two terms are the exchange and 
correlation contributions, respectively. In the correlation term, defined as the difference be- 
tween the total energy and the kinetic and Coulomb contributions, edp) is the correlation 



energy of a uniform liquid of density p obtained from quantum Monte Carlo calculations [14 



The above functional was very successful in predicting the freezing transition in a SL 



WC [T^. It is tempting to extend this theory to investigate multi-layered structures. Ideally, 



one would like to include the correlation energy for the complete structure, which is difficult 



to calculate [0. For typical structures, however, intra-layer correlations are much more 
important than inter-layer correlations [|1^ ; the latter, therefore, can be neglected in a first 
approximation. We write the energy functional of the bilayer system as Eb[p] = Eq[p]+Ej[p], 
where 



Ej[pir)] = jdv' j dv [|r - rf + (f]~"^ [p\v) - po] [p^(r') - po] , (4) 

is the inter-layer Coulomb interaction, and the superscripts A and B label the charge densities 
of the two layers. Ei,[n\ includes (within LDA) the intra-layer correlation, which mainly 
determines the freezing density. In writing Ei, we have neglected the tunneling between 
the layers as well as the finite thickness of the layers, taking the charge density profile in 
the direction normal to the layers as a delta-function in each layer, which is a reasonable 
approximation for d/aB > 1. 

In DFT one usually solves the Kohn-Sham equations self-consistently in order to obtain 
the charge density which minimizes Ef)[p]. A simpler alternative (from the computational 
point of view) is to guess a functional form for p(r) with one or more variational param- 
eters. Therefore, we model the charge density of the electron solid as a sum of gaussian 
distributions centered at the lattice sites R, generated by the primitive vectors ai, 
a2, in one layer, p'^(r) = («/vr) X]Rexp(— a|r — Rp), and at R + c in the opposite layer, 
p^(r) = (a/vr) X]Rexp(— a;|r — R — cp). To minimize the Coulomb energy, c is such that 
each site in one layer sits in the centre of a cell in the opposite layer, a is the localization 
parameter; it is zero in the liquid phase and non-zero for inhomogenous phases, and, for like 
densities, it is the same in both layers. 

We have minimized the energy Eb[p] with respect to the localization parameter a and to 
the lattice structure for fixed density and inter-layer distance [0. Of the several possible 
lattices, those listed in Table I are stable in different ranges of d/aB and n. The resulting 
phase diagram is summarized in Fig. |1|. First, we focus on the freezing transition. For large 
d, i.e., weak inter-layer coupling, the freezing density tends to the single layer result ~ 32. 
In this limit, the BL WC consists of two nearly uncoupled triangular lattices. The two 
lattices are staggered to minimize the residual long-range Coulomb interaction (phase D). 
With decreasing d the formation of an inhomogeneous phase is favoured by a gain in the 
inter- layer interaction and, accordingly, freezing takes place at larger densities, i.e., smaller 
Tg. If d/aB ~ 32, inter-layer interactions become more important and the liquid freezes into 
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a lattice with a larger coordination number relative to sites sitting in opposite layers, namely 
two staggered squares (phase B) or rhombic (phase C) lattices. Along the freezing line of 
phase C, the angle 9 = arccosai • a2 grows continuosly up to 90° with decreasing d. For 
d/aB ^ 9, the phase which freezes from the liquid is phase A. It consists of two staggered, 
rectangular lattices; the aspect ratio (3 = |a2|/|ai| at freezing is larger than 1 and its value 
depends on d. For /? = 1, phases A and B coincide. 

When n is swept through the freezing density, nj, the value of the localization parameter 
a at which the energy minimum is attained leaps from zero (for n > n/) to a finite value (for 
n < Uf) such that ar^ ~ 0.6. This corresponds to a very delocalized charge density so that 
the Lindemann's ratio 7, defined as the root mean square fluctuation around a lattice site 
divided by the lattice parameter, coincides in practice to its upper limit pO[, corresponding 



to a uniform density, which is, e.g., 7 = 0.373 for phase D, and 7 = 0.408 for phase B. 

In Fig. |l]the open circles are the results of Swierkowski et al. 0, who found a singularity 
in the linear response function of the BL electron liquid. These points correspond very well 
with our calculated BL WC freezing boundary. In Ref. |]^ the symmetry of the lattice which 
freezes from the liquid was not identified, but the wavevector at which the singularity occurs 
is q/kp ~ 2.5. Since q/kp = 2.50 for phase B, while q/kp = 2.69 for phase D, it seems that 
the inhomogeneous phase found in Ref. P] is in fact phase B, again in agreement with our 
calculation. 

We next consider transitions between the different solid phases. As seen above, lowering 
n at large d freezes the liquid into phase D. If the density is further decreased, i.e., is 
increased at fixed d, the intra-layer interaction weakens with respect to the inter-layer inter- 
action. This lowers the energy of phase C which finally becomes energetically favoured. Since 
phases C and D belong to different symmetry classes this is a first order phase transition . 
In the region of phase C, the energy minimum is at 7^ 90°, while 6 = 90° corresponds 
to a maximum; on reducing n or d, the energy minimum moves towards 6 = 90°, which 
finally becomes a minimum and the lattice makes a continuos transition to phase B. An- 
other structural phase transition occurs between phases A and B, at sufficiently low density. 



when d/ttB ~ 9. The transition is a continuous one, where the aspect ratio becomes larger 
than 1 at a critical density, which depends on d. An example of the overall evolution of the 
aspect ratio is shown in Fig. ^for d/as = 14; j3 first reaches a maximum and then decreases. 
When rg ^ d/as, however, the present calculation becomes less accurate as inter-layer cor- 
relations beyond the static Coulomb interaction, which are not important near the freezing 
transition and which were neglected here, become more important. Finally, the A/B, 
B/C and C/D transition curves join the freezing curve at three triple points where the solid 
phases and the liquid have the same energy. 

To conclude, we note that, as an alternative to high magnetic fields, one can use particles 
with a small effective Bohr radius, like electrons in silicon layers or holes p2| , to quench 
the kinetic energy. The range of densities for the BL WC in zero magnetic field should be 



within reach in coupled layers of holes grown along high index crystallographic directions [23 



where both high effective masses, to quench the kinetic energy, and high mobilities, to 
prevent defect-induced localization, can be obtained. 
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FIGURES 

FIG. 1. Phase diagram of the electron bilayer. Solid curves: first order phase transitions. 
Dashed curves: continuous phase transitions. Solid dots: triple points. Open dots: freezing 
transition deduced from Ref. [^. 

FIG. 2. Evolution of the aspect ratio /? for phase A at d/aB = 14 obtained by minimization of 
-Efefn] with respect to j3 and a. 
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TABLES 

TABLE I. Lattice parameters of phases A, B, C and D. a is the nearest-neighbour distance. 

For each phase, the primitive vectors ai and a2, tlie inter-lattice displacement c, the reciprocal 
lattice vectors bi and b2, and the charge density Ug are indicated. /3 = |a2|/|ai| and 6 is the angle 
between ai and a2. 

Phase ai/a a2/a c bi/(27r/a) b2/(27r/a) n^a^ 

A (Staggered rectangular) (1,0) (0,/3) (ai + a2)/2 (1,0) (0, 1//3) 1//3 

B (Staggered square) (1,0) (0,1) (ai+a2)/2 (1,0) (0,1) 1 

C (Staggered rhombic) (1,0) (cos 6, sin 6) (ai + a2)/2 (1, - cot 6) (0, sin"^ 6) sin"^ 9 
D (Staggered triangular) (1,0) (1/2,^3/2) (ai+a2)/3 (1,-1/^3) (0,2/^3) 2/^3 



11 



